Below you can find a link to the MATLAB code calculating the advection-diffusion problem with the finite element method with Streamline Upwind Petrov-Galerkin (SUPG) stabilization for the model Eriksson-Johnson problem.

Let us recall first the advection-diffusion problem, the Eriksson-Johnson model problem (described, for example, in [1]), defined on a square domain \( \Omega = [0,1]^2 \) in the following way. We seek \( u \) such as
\( a(u,v)=l(v) \forall v \) where
\( a(u,v) =\int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y) }{\partial x } dxdy + \int_{\Omega} \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ +\epsilon \int_{\Omega} \frac{\partial u(x,y) }{\partial x} \frac{\partial v(x,y)}{\partial x } dxdy +\epsilon \int_{\Omega} \frac{\partial u(x,y)}{\partial y } \frac{\partial v(x,y)}{\partial y } dxdy \)
\( l(v) = \int_{\partial \Omega } f(x,y) v dxdy \)
\( f(x,y)=sin(\pi y)(1-x) \)
is an extension of the Dirichlet boundary condition to the entire area, while \( \beta = (1,0) \) represents the wind blowing from left to right, while \( \epsilon = 10^{-2} \) is the diffusion coefficient. The quantity \( Pe=1/ \epsilon = 100 \) is called the Peklet number, and it defines the sensitivity of the advection-diffusion problem.

In line 848 we setup the diffusion coefficient \( \epsilon = 10 ^{-2} \).
Its inverse is the so-called Peclet number
\( Pe=\frac{1}{\epsilon} \).
SUPG parameters are setup in lines 849 and 850.
Then a vector of nodes is created (always natural numbers counted from zero) \( knot\_x = [0 \quad 0 \quad 0 \quad 1 \quad 2 \quad 2 \quad 2] \) for the trial space,
and the knot vector along x axis from (0,1) interval, \( points\_x = [0 \quad 0.5 \quad 1] \) which will contain the basis of B-spline functions used to approximate the solution to the Eriksson-Johnson problem. The node vectors for the approximating space and the points on which we span the node vectors are defined separately for the x axis and for the y axis.

The code execution is also possible in the free Octave environment.
Download code or see Appendix 5.

The code is activated by opening it in Octave and typing a command
\( advection\_supg\_adapt \)
After a moment of calculation, the code opens an additional window and draws a numerical and exact solution in it.
The code calculates the residuum error
\( Residual norm: L2 0.015808, H1 0.100417 \)
The code also computes the norm L2 and H1 from the solution
\( Error: L2 31.41 procent, H1 38.74 percent \)
and compares to the L2 and H1 standards from the exact solution.
\( Best possible: L2 2.06 procent, H1 28.73 percent \)
You can see that in order to improve the accuracy of the solution, it is necessary to compact the mesh.